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Abstract 

A microscopic model is proposed for the interactions between sickle hemoglobin molecules based on 
information from the protein data bank. A Monte Carlo simulation of a simplified two patch model 
is carried out, with the goal of understanding fiber formation. A gradual transition from monomers 
to one dimensional chains is observed as one varies the density of molecules at fixed temperature, 
somewhat similar to the transition from monomers to polymer fibers in sickle hemoglobin molecules 
in solution. An observed competition between chain formation and crystallization for the model is 
also discussed. The results of the simulation of the equation of state are shown to be in excellent 
agreement with a theory for a model of globular proteins, for the case of two interacting sites. 

1 Introduction 

The condensation of globular proteins from solution is a subject of considerable experimental and theo- 
retical activity. One one hand, it is important to grow high quality protein crystals from solution in order 
to be able to determine protein structure (and thus function) from X-ray crystallography. On the other 
hand, many diseases are known to be caused by undesired condensation of proteins from solution. In 
both cases one needs to have a reasonably detailed model of the protein-protein interactions in solution 
in order to predict the phase diagram, condensation rate and growth kinetics. This is a major challenge 
to theorists, as these protein-protein interactions arise from many sources and are still relatively unknown 
in most cases. 
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An important example of undesired protein condensation occurs with sickle hemoglobin (HbS) molecules 
in solution. It is known that deoxygenated sickle hemoglobin molecules in red blood cells can undergo a 
two step nucleation process that leads to the formation of polymer fibers in the cell jjj H] ■ These fibers 
distort the cells and make it difficult for them to pass through the capillary system on their return to 
the lung. A direct determination of the homogeneous nucleation of HbS fibers in vitro has shown that 
the nucleation rates are of the order of 10^ — 10* cm^^s^^ and that the induction times agree with 
Zcldovich's theory P] . These rates are comparable to those leading to erythrocyte sickling in vivo. They 
are also approximately nine to ten orders of magnitude larger than those known for other protein crystal 
nucleation rates, such as lysozyme. 

Consequently a goal of current research is to understand at a molecular level this nucleation process 
and by controlling the conditions on which the nucleation depends, to slow down the nucleation rate 
such as to prevent the polymerization from occurring while HbS is in its deoxygenated state in the cells. 
To do this requires understanding the protein-protein interactions, in order to predict the phase diagram 
and nucleation rate for sickle hemoglobin molecules. The phase diagram for HbS is only partially known 
experimentally. It is known that there is a solubility line separating monomers and fibers jH Ej and 
evidence exists for a spinodal curve with a lower critical point In a previous publication ^ we 
obtained a phase diagram that was qualitatively similar to this, namely, a liquid- liquid phase separation 
with a lower critical point. In addition, we determined the location of the liquidus and crystallization 
lines for the model, as shown in Figure ^ 

However, although yielding a lower solution critical point, this model was unable to predict the for- 
mation of polymer fibers, as it was based on a spatially isotropic, short range protein-protein interaction 
(e.g. a square well or a modified Lennard- Jones potential energy). Fiber formation clearly requires 
anisotropic interactions. In this paper we propose an anisotropic model for the HbS-HbS interactions, 
based on an analysis of the contacts for HbS crystals from the protein data bank. We also define an 
order parameter to describe the polymerization of this model. As the full model is complex and involves 
several unknown interaction parameters, we study a simplified version of the model (a two patch model) 
in order to gain some insight into the nature of the fiber formation. We determine some aspects of the 
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Figure 1: Phase diagram of a modified Lemiard- Jones model including solvent-solute interactions. Open 
triangles denotes liquid-solidus line; open circles denotes fluid- fluid coexistence. From 7 . Details of the 
model are given in this reference. 

phase diagram for the two patch model via Monte Carlo simulation and biasing techniques and show in 
particular that it yields one dimensional chains that are somewhat similar to the polymer fibers formed 
in HbS nucleation. Real HbS fibers, however, have a diameter of about 21 nm. In addition, the strands 
within the fiber are packed into double strands. Thus the two patch model is too simple to describe 
the polymer fiber phase transition observed in HbS. Future work will be necessary to obtain reasonable 
estimates of the interaction parameters in the full model, in order to obtain a realistic model for the 
polymer fiber phase transition. 

The outline of the paper is as follows. In section 2 we propose an anisotropic interaction model for 
the pair interactions between HbS molecules. In section 3 we define an order parameter that measures 
the degree of polymerization in the system. In section 4 we present the results of our Monte Carlo sim- 
ulation for a two patch approximation to the full model, since we are unable to make realistic estimates 
for the interaction parameters for the full model. A biasing technique is used in order to examine the 
nature of the chain formation. In section 5 we summarize the results of a perturbation theory as applied 
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to an eight patch model and to our two patch model. In the latter case we show that the simulation 
results are is excellent agreement with this theory. In section 6 we present a brief conclusion and suggest 
directions for future research on this subject. 

2 Anisotropic model for the Hemoglobin S polymerization 

Protein molecules in general, and sickle hemoglobin molecules in particular, are very complicated objects, 
typically consisting of thousands of molecules. There are many types of forces between protein molecules 
in solution, such as Coulomb forces, van der Waals forces, hydrophobic interactions, ion dispersion forces 
and hydrogen bonding. Although these interactions are complex, considerable success in predicting the 
phase diagrams of several globular proteins in solution has been accomplished by using rather simple 
models with spatially isotropic interactions. These models share in common a hard core repulsion to- 
gether with a short range attractive interaction (i.e. the range of attraction is small compared to the 
protein diameter). However, in general the protein- protein interactions are anisotropic, often arising from 
interactions between specific amino acid residues on the surfaces of the interacting molecules; i.e., cer- 
tain areas of a protein surface interact with certain areas of another molecule's surface. Thus there 
has been some recent work attempting to model these anisotropic interactions. In such models a given 
protein molecule is represented by a hard sphere with a set of patches on its surface |SI [S] llQl 1111 1121 [T5| . 
Intermolecular attraction is localized on these patches. Typically the models assume that two protein 
molecules interact only when they are within the range of attractive interaction and when the vector 
joining their centers intersects patches on the surface of both molecules. The fluid-fluid diagram for such 
a model was studied by Kern and Frenkel in a Monte Carlo simulation and by Sear theoretically jjl] . 
In these studies all patches are assumed to interact with each other equally. This approximation gives a 
good qualitative picture of the possible importance of anisotropy in protein phase diagrams. However, 
when we consider the behavior of a specific globular protein, this approximation is too simple and does 
not reflect the actual structure of protein aggregates. This is particularly important in the fiber forma- 
tion of HbS molecules in solution. Figure [21 shows interacting HbS molecules, with the different pairs of 
interacting regions in the crystal state of HbS. We use this information, in accordance with the contact 
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Figure 2: Location of the residues participating in different contacts. Yellow areas denote the lateral 
contacts, green areas denote the axial contacts. 

information in the HbS crystal |14[ [T5] , to develop a model for the anisotropic interactions between the 
HbS molecules. 

To develop an anisotropic model to describe these interacting molecules, we allow for the possibility of 
different interaction strengths between different pairs of these interacting patches. To characterize such 
interactions, wc introduce an interaction matrix {eki}mm, where is the strength of interaction between 
the A;*'* and l*^ patches and m is the total number of patches on a protein surface. This is a symmetric 
matrix that describes the strength of interaction between each pair of patches. In our particular model 
we choose a square well attraction between patches, although this is not a necessary restriction on the 
interactions. Thus the interaction matrix consists of the well depth values for the different patch-patch 
interactions. We define the pair potential between two molecules in this case in a way similar to 
but generalizing to the case of different patch-patch interactions; 




(1) 



k,l 
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Here m is the number of patches, fi^ is the orientation (three Euler angles, for example) of the z*'' 
molecule, Ufj is the square well potential for the unit well depth and hn. is the fc*'' patch direction in the 
laboratory reference of frame: 

n,fc = nmnl (2) 

Here n° is the fc*'' patch direction in the i*'* molecule reference of frame and R{Qi) is a rotation matrix 
between the i*'' molecule and the laboratory reference of frame. 8(2:) in is a step function 

{1 if a; > cos(5fc; 
(3) 
if a; < cosdk- 

where 6i is the half open angle of the i*'* patch. In other words, Q{rijhik) is equal to 1 when the vector 
joining two molecules intersects the A:*'' patch on the surface. If patches do not overlap then the sum in 
equation Q has at most one term. The radial square well dependence with range A is given by 

00 for r < (T 
^SW = S -1 forcr<r<Acr 
for r > Act 

Until now what we have done is valid for any pair of molecules interacting through pairwise patch 
interactions. We now consider the specific case of the HbS molecule. As noted above, in agreement with 
contact information for the HbS crystal |14l [TH] , this molecule has two lateral and two axial patches 
that are involved in intra double strand contacts and four more patches involved in inter double strand 
contacts. Thus we have eight possible patches for the HbS molecule (FigureEj). One of the lateral patches 
contains a /?6 valine residue and another lateral patch contains an acceptor pocket for this residue. But 
it is known 'E' that only half the mutated sites are involved in the contacts in the HbS crystal. Thus we 
have another possible set of 8 patches similar to the first one. The total number of patches is therefore 
sixteen (two equal sets of eight patches). The interaction matrix can be built assuming that the first 
lateral patch (of any set) can interact only with the second lateral patch (of any set). The same is true 
for axial patches. The remaining four patches in each set can be divided into pairs in a similar way, in 
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accordance with Jl]. This gives the following interaction matrix (for one set of eight patches): 
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where ei is the strength of the lateral contact, €2 is the strength of the axial contact, £3 and £4 are the 
strength of the inter double strand contacts. We will refer to this model as the "full model" for HbS. 



3 Polymerization order parameter in a system of patchy hard 
spheres 

One of the goals of a study of HbS molecules in solution is to calculate the free energy barrier that 
separates the monomer solution from the aggregate (polymer chains/fibers) state. For this purpose we 
have to specify an "order" parameter that measures the degree of polymerization in the system. The 
structure of the aggregate depends strongly on the configuration of patches. Therefore, to separate 
the aggregate state from the monomer state, the order parameter should reflect the configuration of 
patches. (Note that the order parameter as defined below is only zero in the case in which there are 
only monomers.) Since in our model the regions on the molecular surface not covered by patches do 
not interact (except through the hard core repulsion), we can measure the degree of polymerization by 
measuring the fraction of the patches involved in actual contacts. 

We assume that any two particles at any given time can have no more than one contact between each 
other. This condition is a little stronger than just a non-overlap of the patches. For each pair of particles 
we introduce a quantity that shows how much these particles are involved in polymerization (basically, 



7 



showing the presence of the contact between thcni): 

Np 

il;ij{r,,rj,fli,nj) = ^Wfci/fe(fy ■ hik)fi{-v^j ■ iiji) (5) 

k,l 

where Wki is a weight of the contact between the k*^ patch of the i*'* molecule and the l*^ patch of the j*'' 
molecule. We choose the weight matrix to be the interaction matrix. fk{x) is equal to x for x > cos 5k 
and is zero otherwise. Due to our assumption of only one contact per pair of particles, the sum in Q 
has at most one nonzero term. We next define the order parameter of one particle to be 

Mr^) = (6) 
The term in the denominator is a normalization constant. The order parameter of the whole system is 

^ = ^E^^ (7) 

i 

This choice of order parameter reflects the patch configuration; the magnitude of the order parameter 
increases as the number of contacts in the system increases. It is also rotationally invariant. However, 
this construction has its disadvantages. For some choices of weight matrices it is possible that a fewer 
number of contacts could lead to a larger order parameter if these contacts has significantly larger weights. 
However, the choice of the weight matrix equal to the interaction matrix seems to be natural. 



4 Two-patch model and one-dimensional fiber formation 



The full model described earlier for the sickle hemoglobin molecule is complex and has several interaction 
parameters which have not yet been determined from experiment. Because we are interested in studying 
the fiber formation, we use a simplified model that still can produce fiber chains. We simplify the original 
model by reducing the number of contacts. Since one important feature of HbS fibers is the presence of 
twisted, quasi-one dimensional chains, we consider a system of particles with only two (axial) patches. 
This model is obviously not an accurate description of interacting HbS molecules, but it can lead to the 
formation of one dimensional chains. The interaction matrix for the simplified two-patch model is just 
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^1 Oj 



(8) 




Figure 3: Intermediate state of one dimensional chain formation using the two patch model. The actual 
particle size is equal to the size of the colored spheres. Particles that are not involved in any chains are 
scaled down and shown as small blue circles. Particles that are part of some chains are shown as colored 
spheres. Color shows the depth of the particle (z-axis); blue is the deepest and red is the least deep. 
Simulations were performed at low temperature and high supersaturation. 

Since there are two patches on each sphere, there is only a fluid-solid phase transition. (The full model 
described above can have a fluid-fluid phase transition as well.) The formation of the one dimensional 
chains, therefore, is a gradual transition from the fluid phase as the density of the molecules is increased. 
Figure shows chains that result in our simulation. 

Since the formation of these chains is a gradual transition as one increases the density, they do not 
arise from homogeneous nucleation. A necessary property of nucleation is the existence of a nucleation 
barrier in the free energy dependence on the order parameter. This barrier should separate two wells; 
one well corresponds to the metastable phase, the other to the more stable phase. In our case there will 
not be such a barrier. To examine the nature of the chain formation, we determined the dependence of 
the free energy on the patch order parameter by performing two series of umbrella sampling Monte Carlo 
simulations. The flrst set of simulations was done at T = 0.185 and P = 1.0. 

In this case the initial liquid state does not crystallize in the absence of the biasing and remains in the 
liquid state. The order parameter has only one minimum about "^patch = 5* ^ 0.15, corresponding to a 
mix of monomers and dimers (Figure 0)). As we increase the pressure to P = 1.6 the free energy now has 
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Figure 4: Plot of the free energy Afl in units of kT versus the order parameter for the patches, 4", defined 
in the text. Simulations are at a temperature kT = 0.185 and pressure P — 1.0. The minimum at ^E" 
around .15 corresponds to a liquid state, which is a mixture of monomers and dimers. 

a minimum at a lower value of ^' that corresponds to the liquid state (figure [SJ and a second minimum 
at ^' « 1. This second minimum corresponds to a crystal state with all the patches involved in contacts. 
Thus we see that this double-well free energy describes a liquid-solid phase transition rather than to a 
monomer- (one dimensional) fiber transition. This liquid-solid transition is what one would expect for 
this model. 

Thus, as expected, there is no nucleation mechanism for the monomer-fiber transformation. In the case 
of T = 0.185 and P — 2.35 the system simultaneously crystallizes and increases its number of patch 
contacts. The system successfully reaches an equilibrium crystal state and the free energy has only one 
minimum at this state, as seen in Figure|H| However, at lower temperature the picture is quite different. 
At r = 0.1 and P — 0.01 the fibers form before crystallization can occur. The free energy has one 
minimum at 5" around 0.9, but the system is not crystallized. As the set of fibers is formed, the dynamics 
slows down significantly and the system becomes stuck in a non-equilibrium state. Figure |21 shows an 
example of a typical configuration for this noncquilibrium state, corresponding to a set of rods in a 
"glassy" state. 

The umbrella sampling simulations were performed on a system of N particles, with = 500, in the 
NPT ensemble, with a range of interaction given by A = 1.25. The equation of state simulations were 
also performed in the NPT ensemble with 500 particles. The details of the umbrella sampling technique 
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Figure 5: Plot of the free energy Afi in units of kT versus the order parameter for the patches, Vf, 
defined in the text. Simulations are at a temperature kT = 0.185 and pressure P = 1.6. The minimum 
at ^' around 0.2 corresponds to a mix of monomers and dimcrs, while the minimum at ^! close to 1 
corresponds to a crystal state in which all the patches are in contact. 
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Figure 6: Plot of the free energy Aft in units of kT versus the order parameter for the patches, 4*, defined 
in the text. Simulations are at a temperature kT — .185 and pressure P = 2.35. The minimum at ^' in 
the vicinity of 0.9 is the crystal state. 
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can be found in |16[ll7f . In short, this method is based on a biasing of the original interactive potential 
in such a way that the system is forced to reach otherwise inaccessible regions of the configuration space. 
In particular, for simulations starting in a liquid state, the system is forced towards larger values of the 
order parameter. Since the actual dependence of the free energy on the order parameter (for the entire 
range of values of 5") is not known, the biasing function is chosen to be quadratic: 

where the parameters k and determine which region of values of ^' would be sampled in the simulation. 
By changing these parameters we can sample the entire region of values of An interesting result of 
the simulations is that at intermediate pressures the system when biased to large values of 5* starting 
from an initial liquid state has a very large volume. However, if the system is started from an initial (fee) 
crystal state, the volume of the system remains small (still in a crystal state), while the order parameter 
value is around 1. This observation suggests that at not very high supersaturation the biased system 
started from the fluid ends up by forming a few long flbers, rather than a set of flbers that are packed 
into a crystal lattice. The fiber formation time therefore is much smaller than the crystallization time 
at low supersaturation. Starting from aa fee initial condition, however, the particles just reorient within 
the crystal lattice to form the fibers, remaining in the crystal state. At higher pressure, and therefore 
higher supersaturation, (such as in Figure [SJ, the crystallization occurs in a time comparable with the 
fiber formation time. While the system is forced to form the fibers, these fibers pack into the crystal 
lattice. 

4.1 Equation of State 

The equation of state for the two patch model is shown for the two low temperatures studied in Figures 
□ and El 

Also shown in these figures are results from a theory for the m site model ® of globular proteins due 
to Sear ^T]. This model is identical to our patch model defined in section 2 for the case in which the 
various interaction parameters are equal. Sear studied the m-site model (m patches) using the Wertheim 
perturbation theory ^ for the fluid phase and a cell model for the solid phase ^H] and showed that 
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Figure 7: A plot of the pressure, p, versus density, p, in the fluid phase at kT = .17( in units of the well 
depth e). The open circles are the results of the simulation. The asterisks are theoretical results obtained 
from Sear [III. The number of patches m = 2, with a patch angle of about 6 — 52 degrees and a range 
of interactionA = 1.25, i.e. Vf. — 1.25cr, where a is the hard core diameter. 
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Figure 8: A plot of the pressure versus density in the fluid phase at kT — .185 The open circles are the 
results of the simulation. The asterisks are theoretical results obtained from Sear Same values for 
parameters as in Figure [7| 
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Figure 9: Theoretical prediction for the phase diagram for the two patch model discussed in the text. 
The theory is due to Sear JI]- Same values for parameters as in Figure El 

the model exhibits a fluid-fluid transition (for m > 2) which is metastable with respect to the fluid- 
solid transition for most values of the model parameters. For m — 2, however, there is only a fluid-solid 
transition. As can be seen from Figures|7|and|Sl the theory yields results which are in excellent agreement 
with the results of our simulation. For completeness we show the theoretical prediction for the fluid-solid 
transition for m = 2 in Figure |51 assuming a fee crystal structure. 

We also show in Figure EH the results of the theory for the case m = 8, as the model discussed in 
section 2 has 8 pairs of interacting patches. Since the interactions between the different sites in the model 
studied by Sear are assumed to be equal, the model lacks the anisotropy discussed in section 2 that is 
necessary to account for the fiber formation in HbS molecules. Nevertheless, it is quite instructive to 
know the phase diagram for this case. 

The fluid-fluid binodal curve has an upper critical point for m > 2 (e.g. Figure unlike the case 
for HbS. In that case experimental measurements by Palma et al ^ display a spinodal curve. Such a 
spinodal implies the existence of a binodal curve with a lower critical point . However, as shown by 
Shiryayev et al the lower critical point reflects the crucial role of the solvent in the case of HbS in 
solution. The solvent is not taken into account in our model, but presumably if one would include a 
solute-solvent coupling similar to that of 7 , this coupling could change the phase diagram shown in 
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Figure 10: Theoretical prediction for phase diagram for the model of HbS discussed in the text, in which 
all interaction parameters are equal. Theory is due to Sear Here m — 8, X — 1.05 and 5 — 51 

degrees. The fluid-fluid transition is metastable. 

Figure^lto one with a lower critical point, as found by Shiryayev et al 0, e.g. Figure^ 

Finally, we note that Jackson et al [HI used Wertheim's theory for the two site model to predict that 
the fraction of molecules that are present in chains of length n is given by 



nX^{l-X) 



n-l 



(9) 



while the average change length, jn^, is given by 

< 72 >= 1/X. 

Here X is the fraction of sites that are not bonded to another site and is given by 

X 



(10) 



(11) 



l + [l + 4pX5£,exp(/3e)]i/2 

where g^j, is the contact value of the pair distribution function for a fluid of hard spheres. The quantity 
K is given by the expression 

ii: = 7rCT2(r^_cr)(l-cos(2(S))2. (12) 

For example, the theoretical prediction for the fraction of dimers and the average chain length as a 
function of density (at kT — .185) is shown in Figures ITTland lT^ respectively. Also shown for comparison 
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Figure 11: Theoretical prediction for the fraction of dimers as a function of density at kT — .185. From 

m 

in Figure 1121 is an approximation for the average chain length obtained from our simulation results for 
the order parameter. Had we chosen to use a function / = 1 in our definition in section 3 (rather than 
f{x) — X for X > cos (5) the order parameter would have been equal to 1 — X, i.e. the fraction of sites 
that are bonded to another site. In that case the average chain length would be equal to 1/(1 — ^P). As 
shown in the figure, 1/(1 — ^) is in general less than the average chain length, due to our choice of f in 
section 2. 



5 Conclusion 

The results obtained in the previous section raise several important questions. For some thermodynamic 
conditions the system crystallizes and the fibers align along each other to form a fee like structure. 
For other thermodynamic conditions the fiber formation prevents the system from crystallizing and it 
remains in a non-equilibrium glassy state. Is there some boundary between these two behaviors? In 
protein crystallization science a somewhat similar phenomenon is known as " gelation" . The boundary 
between successful crystallization and the nonequilibrium "gel" state is the "gelation line". (This line 
is obviously not an equilibrium phase boundary.) This "gel" state is not a gel in the usual sense. It is 
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Figure 12: Theoretical prediction for the average chain length as a function of density at kT = .185 
(open circles). Also shown are the simulation results for 1/(1 — 5') (asterisks)-see discussion in text. 

more like a glassy state in which the aggregates form, but due to their formation the dynamics slows 
down significantly and the aggregates cannot subsequently form crystals. Thus the system becomes 
stuck in this glassy state. Recently mode coupling theory has been used to predict the gel line for protein 
crystallization ;19 . It is possible that one can use the same approach for the two patch model. 

Another interesting question arises when we compare the behavior of the two-patch model and real 
sickle hemoglobin . Ferrone's studies WJD show that the formation of the fourteen strand fibers occurs via 
a nucleation mechanism. Molecules aggregate into an initially isotropic droplet which subsequently be- 
comes an anisotropic fiber. The kinetics of the transformation of the isotropic droplet into an anisotropic 
fiber is not well understood. It is believed that this happens through the attachment of molecules to 
active sites of the molecules in the fiber. At some fiber diameter the number of active sites on the fiber 
surface is not sufficient to induce layer by layer growth of the fiber in the direction perpendicular to the 
fiber axis. Thus molecules that attach to the active sites form a droplet that then detaches from the 
fiber, with the subsequent formation of a new fiber (i.e. the new fiber forms via heterogeneous nucleation 
from the original fiber). This process explains why an anisotropic fiber does not continue to increase its 
diameter. 

The point of this discussion is to note that one of the main differences between the two-patch model 
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and sickle hemogloblin molecules is that the latter forms fibers through nucleation, whereas the former 
does not. One way to improve the two patch model is to increase the number of patches. If one includes 
an additional two active patches, corresponding to the lateral contacts, then most likely this system 
would form non-interacting double strands. If so, this would not lead to any qualitative difference with 
the two patch model. Another approach would be to add several relatively weak patches around the 
particle which can represent inter strand interactions. However, this would not necessarily give two 
distinct nucleation mechanisms (monomer-fiber and fiber-crystal). Rather, it is more likely that the 
crystallization would be anisotropic. In order to produce a nucleation from monomers to fibers it might 
be necessary to have a particular distribution of patches such that at some radius of the anisotropic 
droplet (pre-fiber) its growth in the radial direction is significantly depressed. 
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